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f») . Abstract 

In a previous work [arXiv: 1009.4363], we have studied the evolution of a scalar field with 

r^ i ' a quartic coupling, driven by a classical source that initializes it to a non-perturbatively large 

value. At leading order in the coupling, the evolution of this system is given by classical 

£SJ ■ solutions of the field equation of motion. However, this system is subject to a parametric 

^ ' resonance that leads to secular divergences in higher order corrections to physical observables. 

We have proposed a scheme that resums all the leading secular terms: this resummation leads 

to finite results at all times, and we have observed also that it makes the pressure tensor of 

the system relax to its equilibrium value. 

In the present paper, we continue the study of this system by looking at finer details of its 
£ ■ dynamics. We first compute its spectral function at various stages of the evolution, and we 

observe that after a fairly short transient time there are well defined massive quasi-particles. 
We then consider the time evolution of the momentum distribution of these quasi-particles, and 
we show that after a stage dominated by the parametric resonance, this distribution slowly 
evolves to an equilibrium distribution. Interestingly, this distribution develops a transient 
chemical potential, signalling the fact that number changing processes are much slower than 
the elastic ones. 



1 Introduction 

The issue of thermalization is a very challenging question in the theory of heavy ion collisions. On 
the one hand, a lot of circumstantial evidence suggests that the matter formed in these collisions 
behaves like a nearly perfect fluid at very short times after the collision (this is based on a com- 
parison of flow measurements with predictions of hydrodynamical models - see [1] for instance). 
This fact is commonly interpreted as a sign of the fast thermalization of the quark-gluon matter 
produced in nucleus- nucleus collisions 1 . On the other hand, justifying this thermalization from 



1 Note however than some systems with disordered field configurations may also exhibit an anomalously low 
viscosity [2, 3]. 



first principles in QCD has proven to be very difficult, and to this date the question has not been 
satisfactorily answered. 

The main issue is the following: a semi-hard momentum scale -the saturation momentum Q s - 
is dynamically generated by the non-linear gluon interactions that occur in nuclear collisions at 
high energy [4, 5, 6]. This scale controls the so-called gluon saturation phenomenon, that alters 
the gluonic content of a nucleus at high energy, but it is also this scale that controls the final 
state interactions that occur after the collision, at least initially. However, the saturation scale 
Q s increases with the energy of the collision, and therefore the relevant QCD coupling effectively 
decreases. It thus seems that the relevant degrees of freedom are rather weakly coupled in collisions 
at very high energies, rendering difficult the obtention of short thermalization times [7]. 

To compute the initial gluon production in heavy ion collisions at high energy, we have at our 
disposal an effective theory -the Color Glass Condensate [8, 9, 10, 11, 12, 13, 14]- that consistently 
takes into account the non-linear saturation effects. In this effective theory, the degrees of freedom 
that populate the light-cone wave-function of a nucleus are divided in two classes, according to 
their longitudinal momentum. Partons with a longitudinal momentum above a certain cutoff A 
are approximated by a collection of static (because they appear highly boosted in the frame of the 
slower partons) color sources of density p a {x±), while the partons with longitudinal momentum 
below A are treated as standard gauge fields 2 A^(x). In this effective theory, the sources p a are 
stochastic variables (reflecting the fact that we do not know the precise arrangement of the fast 
partons at the time of the collision), with a probability distribution W[p]. 

From the above description of the CGC, it may seem that observables could depend on the 
unphysical cutoff A, that has been introduced by hand as a separation scale between the fast and the 
slow partons. This worry becomes manifest when one computes loop corrections in this framework: 
they usually contain logarithms of the cutoff A. Fortunately, one can show that for sufficiently 
inclusive observables in the collision of two nuclei, these logarithms are universal in the following 
sense: (i) they do not depend on the observable under consideration, but are instead properties of 
the initial state, and (ii) each logarithm of A can be assigned to one of the two projectiles (meaning 
that its coefficient does not depend on the nature of the second projectile). This universality allows 
to absorb all these logarithms into a redefinition of the distributions W[p] of the two projectiles. 
After this is done, these distributions depend on the cutoff A (in such a way as to cancel the A 
dependence that arises from loop corrections), and their changes under variations of the cutoff are 
controlled by the JIMWLK evolution equation [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. 

The Color Glass Condensate, supplemented with the JIMWLK equation, therefore provides 
a consistent framework for computing observables relevant to heavy ion collisions, such as the 
spectrum [26, 27, 28, 29, 30, 31, 32] of produced gluons or components of the energy-momentum 
tensor [33]. However, these computations have also shown that the gluonic matter produced in 
these collisions is rather far from local equilibrium [34]. In fact, its pressure tensor is highly 
anisotropic, with two positive transverse pressures and a negative longitudinal pressure. Such a 
form of the pressure tensor is arguably too far from the equilibrium form to justify the applicability 
of hydrodynamics (even with viscous corrections). 

This computation, a leading order calculation in a s improved by the resummation of all the 
leading logarithmic contributions discussed above, is however not the end of the story. It has been 
noted in various situations [35, 36, 37, 38, 39, 40, 41, 42] that the classical solutions of Yang- 



2 At the leading order of this description, the gluons are dominant over the quarks. Only the former are taken 
into account in the description. 



Mills equations that constitute the LO answer in the CGC approach are unstable under small 
perturbations. Practically, this means that loop corrections contain secular divergences -i.e. terms 
that have an unbounded growth with time- and that the perturbative expansion breaks down 
after a certain time (this time can be estimated to be of the order of Qj 1 , up to logarithms). An 
additional rcsummation is therefore necessary in order to circumvent this problem and to obtain 
expressions that remain well behaved at all times. 

Interestingly, the structure of the NLO (1-loop) corrections to inclusive observables in the CGC 
formalism suggests a natural resummation scheme, that in a loose sense amounts to exponentiating 
the 1-loop result. Such a resummation first appeared in the context of inflation in cosmology 
[43, 44]. In the context of the color glass condensate and heavy ion collisions, it was sketched in 
[45] in the case of gauge fields, and more thoroughly justified 3 in the case of scalar fields in [48] . The 
reason why a </> 4 scalar field theory in four dimensions is an interesting playground for testing these 
ideas is that this theory, like Yang-Mills theory at the classical level, is scale invariant, and that 
its classical field configurations are also subject to an instability (due to parametric resonance 4 ). 
The rcsummation wc advocated in [48] can be formulated in two equivalent ways: its derivation 
leads first to a rather formal expression, in which the rcsummation is expressed as the action of 
the exponential of a certain operator on the LO result. This form is not practical for a numerical 
implementation, but is equivalent to averaging the LO order result (which is expressed in terms 
of a classical field) over a Gaussian ensemble of initial conditions for the classical field. We then 
showed that the proposed resummation indeed leads to observables that remain finite at all times, 
therefore solving the original problem of secular divergences. Moreover, by studying the time 
evolution of the pressure tensor, we showed that after the resummation, the pressure converges 
towards its equilibrium value in a fairly short time -something that would not happen with the 
unresummed expression. 

In this paper, we pursue the study of this system of scalar fields, in order to elucidate the 
microscopic state of the system. In particular, we would like to know whether the system is 
thcrmalizcd when its pressure has relaxed to its equilibrium value, or whether on the contrary one 
could have an equilibrium- like pressure tensor while the system is still far from equilibrium 5 . A 
natural quantity to study in order to address this question is the occupation number in the system, 
and its time evolution. However, even before computing the occupation number, it is interesting to 
ask whether the system can be described in terms of quasi-particles (this is not trivial: although the 
system is weakly coupled, it is also very dense, and strong collective effects may render the quasi- 
particles completely unstable) . We start with a reminder of the model and of the main results of [48] 
in the section 2; then we compute in the section 3 the spectral density of the system, after having 
justified that it can be resummed in the same way as the energy- momentum tensor. In the section 
4, we continue our study with the occupation number. We compute it as a function of momentum 
up to large times, and identify several stages in its time evolution. In particular, we discuss the 
possible occurrence of turbulence and Bose-Einstein condensation. From the occupation number, 
we perform in the section 5 several tests of the quasi-particlc picture (e.g. compare the measured 



3 This rcsummation was also derived in a completely different approach in [46]. The full spectrum of initial 
fluctuations in Yang-Mills theory in the Fock-Schwinger gauge has been recently derived in [47]. 

4 The problem of secular divergences and parametric resonance has been addressed in a number of other approaches 
in Quantum Field Theory, including perturbation theory (in the regime where the resonant modes arc still small) 
and the 2-Particlc Irreducible rcsummation [49, 50, 51, 52, 53, 54]. 

5 In the 2PI approach to the relaxation of a linear sigma model, it has already been observed that an equation of 
state can be obtained much earlier than the complete thermalization of the system [55]. 



mass of the quasi-particles, with its value at 1-loop, including medium effects) and compute how 
the entropy of the system evolves with time. A summary and further discussion can be found in 
the section 6. In the appendix A, we discuss some aspects of the classical dynamical system to 
which the quantum field theory is equivalent in our resummation scheme. Some basic properties 
of the Liouville equation are recalled in the appendix B. 

2 Setup of the problem - Summary of previous results 

2.1 Model and resummation of secular terms 

In [48], we considered a real scalar field theory, with quartic self-interactions, and coupled to an 
external source. The Lagrangian of this model reads 

£^(^)(c^)-^ 4 +J<£. (1) 

In this model, the source J is nonzero only at negative times, and its purpose is to initialize the 
fields to a configuration that has a large expectation value at t = 0. At positive times, J = and 
the fields evolve by themselves. 

In [48] , we have studied the time evolution of the energy-momentum tensor of the system. At 
leading order (tree level) in the coupling g 2 , it is simply the energy-momentum tensor of a classical 
field: 

T%{x) = d^d'tp-g^^d^) 2 -^ 

□^ + |]V = J, lim <p(x°,x)=0. (2) 

O ! x" — ¥ — oo 

However, at NLO (1 loop) we observed that T^ v is plagued by secular divergences that originate 
from an instability of the above classical solution, due to parametric resonance. In order to cure 
this problem, we developed a resummation scheme that collects the leading secular terms at each 
order of the expansion in g 2 . Our resummed expression reads 

TZn^) = ^v[jd 3 ufi- T u + l -J d 3 ud 3 v J '^0^[ a+k ■ T u ][a„ k ■ T v ]]t£(x) . (3) 

In this formula, T u is the generator of shifts of the classical field on the x° = hypersurface, 
defined formally by 

A A 

a-T u = a(0,u)— — _ + d(0,«)— -j- -r , (4) 

where (fio(u) and <fio(u) are the values of the field and its first time derivative on the surface x° = 0. 
This shift operator possesses a very useful property: if a(x) is a small perturbation on top of a 
classical field ip(x), one has 

a(x) = / d 3 u [a ■ T u ] tp{x) . (5) 



The fields a± k are small perturbations propagating on top of the classical field cp, that are plane 
waves at x° — > — oo, and (3 is the 1-loop correction to ip, 



a + v"(cp) 

U + V"{<p) /3 = 



a±k = 
1, 



lim a± k (x) = e ±ik ' x , 
d 3 k 



{2n) 3 2k 



lim /3(x) = 



X u —¥ — OO 



(0) 



(See [48] for more details.) A crucial result is that eq. (3) is equivalent to a functional integration 
over Gaussian fluctuations of the classical field at a; = 0, 



TZim = J[Da(x)Da(x)} F[a, a] T%[cp + f3 + a] , 



(7) 



where T££ [cpo + f3 + a] denotes the LO energy-momentum tensor evaluated with a classical field 
whose initial condition at s° = is ipo + f3 + a (and likewise for the first time derivative). The 
distribution F[a, a] is Gaussian in a(x) and a(x), with 2-point correlations given by 



(a(x)a(y)) = 
(a(x)a(y)) = 



d 3 k 



a +k {0,x)a- k (0,y) , 



(2ir) 3 2k 

d 3 k 

—3— d +fc (0,a;)d_ fe (0,y) . 



(8) 



(2n) 3 2k 
We refer the reader to [48] for more details about the model and its practical lattice implementation 

2.2 Relaxation of the pressure 




Figure 1: Relaxation of the pressure tensor towards the equilibrium value. All the numerical results 
presented in this paper have been obtained with a coupling constant .9=1. 

The main result of [48] is that the fluctuations resummed in the formula (7), in addition to 
making the resummed quantities free of secular divergences, cause the pressure to relax towards 



its equilibrium value (one third of the energy density for a </> 4 theory in four dimensions). Without 
these fluctuations, the pressure is simply that of a classical field configuration, and it oscillates 
periodically around the equilibrium value, but never relaxes: therefore there is no one-to-one 
relationship between the pressure and the energy density, and thus no equation of state at LO. 
The fluctuations of the initial condition for the classical field in eq. (7) produce a decoherence that 
dampens the oscillations of the pressure, as illustrated in the figure 1, and after a finite time the 
pressure converges towards the value required by the equilibrium equation of state (3p = e). 

In [48] , we have also studied the fluctuations of the energy in a small subsystem. Of course, the 
total energy in the system is conserved and therefore does not fluctuate, but its value in a small 
subvolume can vary due to exchanges with the surroundings. We have observed a rapid transition 
from rather narrow Gaussian fluctuations at the starting time to broader fluctuations that are 
much closer in shape to canonical equilibrium fluctuations at late times. This transition appears 
to occurs approximately at the same time as the relaxation of the pressure. 

2.3 Lattice setup for the study of thermalization 

However, the study performed in [48] does not tell us much about the microscopic evolution of the 
system: indeed, since the energy density and the pressure of the system are intensive quantities 
that sum over all the modes of the system, the existence of an equation of state identical to the 
equilibrium one does not imply that the system is in full equilibrium. Equilibrium requires a 
much more stringent microscopic arrangement of the system, in which the energy is distributed 
among the various modes in a very specific way. In the present paper, we wish to pursue this 
study by looking at more microscopic aspects of the system: namely the existence of quasi-particlc 
excitations, and the evolution in time of their momentum distribution. In particular, we would 
like to know how the energy of the system, initially introduced in the zero mode via a spatially 
homogeneous external source, eventually goes into higher momentum modes. 

A crucial ingredient in this process is the parametric resonance that exists in the </> 4 scalar field 
theory, and therefore it is important that the ultraviolet lattice cutoff be large enough to comprise 
the resonance band. If the source J is parametrically 

O 3 

J ~T' (9) 

then the resonance band is located at momenta of order k ~ Q, and the UV cutoff A must therefore 
obey Q < A. Setting up the lattice cutoff in this way is sufficient to study the evolution of the 
system at short times, because on these time scales the occupation number remains small above 
the resonance region, as we shall see later in this paper. 

This is however not sufficient if we want to study the approach of the system to thermal 
equilibrium. In order to see this, recall that the energy density in the system is parametrically 

O 4 

If thermal equilibrium is achieved, this energy density must also be given by the Stefan-Boltzmann 
formula (at least for reasonably weak couplings) , 

e~T\ (II) 



which tells us that the system would equilibrate at a temperature 

T~-^-. (12) 

V9 

For a numerical simulation to be able to approach the equilibrium state, the lattice ultraviolet cutoff 
must be large enough to include modes of the order of the temperature, which is a more stringent 
constraint than simply having the resonance band below the cutoff. At weak coupling, this implies 
that the resonance band should be located towards the soft sector of the lattice spectrum, i.e. in 
a region where the lattice mode density is rather low. In order to still have a significant number 
of lattice modes inside the resonance band, we used a larger lattice (of size 20 3 for the results 
presented in this paper) , and we have chosen the value of the parameter Q so that the resonance 
band is located near fcssl (in lattice units, where the ultraviolet cutoff is at A = ^/\2). 

This choice of Q is significantly lower than the value used in [48] (where we were only interested 
in the early stages of the time evolution, dominated by the resonant modes) . This means a lower 
energy density, and larger time scales. Indeed, since our system is scale invariant, energy density 
scales like Q 4 and all the times scale like 1/Q- The resulting time evolution of the pressure, for 
this choice of Q and a coupling constant g = 1, is shown in the figure 1. 

3 Spectral function and quasi-particles 

3.1 Definition and leading order 

Before we study the time evolution of the occupation number in the system, it is interesting to 
ask an even more elementary question: can the system be described in terms of quasi-particles, 
or on the contrary does it interact so strongly that no identifiable quasi-particles show up in its 
spectrum? To that effect, one can compute the spectral function, defined as the imaginary part of 
the Fourier transform of the retarded propagator 6 : 

p(oj,k;y°) = 2Im f dtd 3 x e tujt e - tkx G R {y° + t,x,y° 7 0) . (13) 

b 

In this formula, the retarded propagator is normalized so that 

□, + V"bp{x)j\ G H (x, y) = S(x - y) (14) 

for a classical field configuration ip. A system in equilibrium is invariant under translations in 
time, and therefore its spectral function defined in this way is in fact independent of the time y°. 
However, for transient systems that are not yet in equilibrium, the spectral function will evolve 
with time and the y° dependence is important. 

At leading order, we simply need to obtain the retarded propagator in a classical background 



6 We assume here that the system is spatially homogeneous. This is the case in our setup since the source J does 
not depend on x. 



field tp(x), 



G\°{x,y) = Y, 




(15) 



where the grey blobs denote the retarded classical field f(x) and the lines with an arrow the bare 
retarded propagator (the sum is over the number of insertions of the classical field, from to +co). 
In a numerical calculation, the simplest way to compute this propagator is to consider a small field 
perturbation a(x) in that background, that obeys the following equation of motion 

a x + V"((p(x)j\ a(x) = . (16) 

The fluctuation a(x) can be related to its value at the time y° by the following Green's formula 



a(x) 



d 3 y 



G%>(x,y) {d° y a{y)) - (tfG%>(x,yj) a(y) 



(17) 



that involves precisely the propagator we are looking for. Thus, by choosing the following initial 
conditions at time y , 

a(y°,y) = 0, d° y a(y°,y) = 5(y) , (18) 

the perturbation a(x) is precisely the propagator we need in eq. (13), 



G^(x°,x,y°,0) = a(x t x) 



(19) 



Thanks to this observation, we reduce the problem of finding the retarded propagator at leading 
order to that of solving the equation (16) with the initial conditions of eq. (18), which is easily 
doable numerically 7 . 



3.2 Next to leading order and resummation 

At next to leading order, three topologies must be evaluated 8 : 






(20) 



where the both the propagators and the vertices are dressed by the classical field f{x). The first 
topology is the same as the retarded propagator at leading order, in which one of the (p insertions 
has been replaced by a 1-loop tadpole j3 defined in eq. (6). Given that this tadpole is given by (see 
the eq. (39) in [23]) 



P(x) = 



d 3 u (3 ■ T z 



- I d 3 ud 3 v 
2 



d 3 k 



(2n) 3 2k 



[a+k -Tulla-k -T„] <p(x) , 



(21) 



7 On a lattice, the delta function that appears in the initial condition for c9jJa(j/°, y) becomes a Kronecker symbol: 
the derivative is zero at all points of the lattice except at the origin (0,0,0) where it is equal to one. 

8 Although these diagrams seem to involve cubic vertices, this is not the case. These vertices where three lines 
merge are in fact proportional to V"'((p(x)), and are thus proportional to an extra <fi(x) that does not appear 
explicitly in the diagrammatic representation. 



it is easy to check that this contribution is related to the leading order one by the following 
functional identity 



If f d 3 k 

GNLoi^y) = J d 3 uf3 . Tu + _ d ;i ud 3 v / 2fc [ a+fc .T u ][a- fc -T„] 



G£(x,y), (22) 



where the subscript 'same ip' indicates that the two operators T U T„ in the second term should 
act on the same field ip, i.e. 



T U T V 



J same ip 

By using the formula 



tp(xi) ■ ■ ■ <p(x n ) = ^2 ^( x i) ' ' ' ffa-i) 



T U T V 



f(Xi) 



G\°_(x,y) 



d 3 k 



(2w) 3 2k 
the second topology can be written as follows, 

d 3 k 



a+k(x)a^ k (y) , 



ip(x i+1 ) ■ ■ ■ cp(x n ) . (23) 



(24) 



GP™*(x,y) 



d 3 ud 3 v 



(2ir) 3 2k 



[a +k ■ T u ][a_ k ■ T„] 



same V"{<-p) 



G™(x,y), (25) 



where the subscript 'same V"(ip)' indicates that the two operators T^T,, should act on the same 
compound V"(<p) -one operator on each field of V"{tp). The third topology can first be written as 



(26) 



G% h ° 3 (x,y) = J d 4 wd 4 zG™(x,w)^°° p (w,z)G%°(z )V ) , 
where £ 11oo p is the 1-loop retarded self-energy, 



Next, one can rewrite this self-energy as 



-illoop 



(w,z) 



illoop 



(w,z) 



1. 



5^ loop (w,z) = -V'"(<p{w))V'"(<p(z))G i R {w,z) G 1 f_(w,z) + G J i° + (w,z) 



(27) 



(28) 



where the prefactor 1/2 is the symmetry factor of the loop. By combining cqs. (26) and (28) and 
by using (24), one can finally prove 



g-: 



l (x,y) = 



rl 



d 3 ud 3 v 



d 3 k 

{2n) 3 2k 



[a+k ■ Tr„][a_fc • T v ] 



distinct ip' s 



Gf{x,y), (29) 



where the qualifier 'distinct cp's' indicates that the two operators T„T B must act on two fields ip's 
that are inserted at different points on the LO propagator G^° . Adding eqs. (22), (25) and (29) 
therefore simply lifts any restriction on the action of these operators, and we obtain 



G^°(x,y) 



d 3 u (3 ■ T u + i f d 3 ud 3 v J Jy 2k [a+k ■ T„] [a_ fc • T„]] G\° (x, y) . (30) 



This formula is formally identical to the formula we have obtained previously for the energy- 
momentum tensor at NLO, and it leads to the same pathologies due to the presence of secular 



divergences. Likewise, the problem can be cured here by performing the same resummation as in 
the case of the energy-momentum tensor, that amounts to exponentiating the quadratic part of 
the operator in the square brackets in eq. (30) (as in eq. (3)): 



GT UB 



l (x,y) = cxp 



rl 



d d ud 3 v 



d 3 k 



(2ir) 3 2k 



[a +k -T u }[a- k -T v ] G%°(x,y) 



(31) 



This resummation amounts to a functional average over Gaussian fluctuations of the initial condi- 
tion of the classical field at x° = 0, 



grssummed = / [Da(x)Da(x)] F[a,a] G™ [cp + a] 



(32) 



where the Gaussian distribution F[a, d] is defined in eq. (8). Therefore, in order to compute the 
resummed retarded propagator, we should repeat the procedure outlined in the section 3.1 for 
every classical field ip obtained from an ensemble of initial conditions (p + a, where a samples the 
Gaussian distribution ,F[a, d]. 

3.3 Numerical results 



p(o),k) 




Figure 2: Spectral function p(u, k; y° = 0.0) at the initial time. The computation is done on 
a 20 3 lattice, for a coupling constant g = 1. In this plot, k denotes the lattice momentum, i.e. 



cos(2ttI/L) 
in the range [0, L — 



- cos(2irm/L) — cos(2irn/L)) on a L 3 lattice (l,m,n is the triplet of integers 
1] that labels a given momentum state). 



At the initial time (see the figure 2), the spectral function has a fairly complicated structure. 
Although the large k modes have a single spectral peak atw« |fe|, the situation is richer in the 
soft sector. There, besides the main branch that continues to large k, the spectral density exhibits 



10 



additional branches. One of them corresponds to a higher mass excitation, and another one has 
a mass comparable to the main branch but an anomalous dispersion such that the frequency 
decreases while the momentum increases. Therefore, at early times, the quasi-particle picture is 
not a good description of the degrees of freedom in the system. 





Figure 3: Spectral function p(uj,k;y°) at the times y° = 400 (left) and y° = 3000 (right). The 
computation is done on a 20 3 lattice, for a coupling constant g = 1. 

As the time increases, these extra branches in the spectral function decrease in amplitude 
and eventually disappear, starting with the higher mass excitation. Ultimately, only the main 
excitation remains, as one can see in the plot on the right of the figure 3 at a time y° = 3000. At 
intermediate times (such as y° = 400, represented on the plot on the left of the figure 3), one gets 
closer to the spectral function of a system made of quasi-particles, with only small remnants of 
the structures that existed at early times. It is interesting to note that the characteristic time for 
the disappearance of the extra branches in p(w, k; y°) is comparable to the relaxation time of the 
pressure, that we have found in the previous section to start at a time of the order of y° ~ 100. 

3.4 Quasi-particle mass 

In order to further assess the existence of quasi-particles in the system, one can try to fit the main 
branch of the spectral function by a function of the form w = yk + to 2 . The result of this fit is 
shown in the figure 4. One sees that the mass to resulting from this fit is not stable until a time 
y a ss 100, and becomes much more regular afterwards. This is in agreement with the previous 
qualitative observation that only the main branch of the spectral function survives after this time. 
Moreover, after y° > 1000, the mass of the quasi-particles that populate the system decreases 
slowly with time, indicating that the system is not yet completely equilibrated (the change in the 
mass of the quasi-particles reflects a change in the occupation number of the various modes of the 
system, that we will study more directly in the following section). If one takes as a crude estimate 
the Hard Thermal Loop [56, 57] expression of the medium-generated mass (see for instance [58], 
PP 41-45), 

d 3 k 

(27rp2/c 

as a function of the occupation number /&, we can interpret the decrease of the mass as a shift of the 
occupation number from low k to higher fc's if the total number of quasiparticles (TV ~ J d 3 k /*.) 

11 



i 2 - 

HTL 



™~~=9 2 T^^r/fe. (33) 




Figure 4: Green line: quasi-particle mass obtained by a fit of the main dispersion branch with a 
function of the form lo = yk + m 2 . Blue line: 1-loop analytic calculation from the occupation 
number. Red line: 1-loop gap equation that resums recursively all the daisy diagrams. (See the 
text in section 5 for explanations regarding the curves labelled m 2 



,2 ,tl and m i-ioo P -; 



is approximately constant. 

A word of caution should be added about the width of the quasi-particlcs. The width of the 
spectral peak in the figures 2 and 3 is probably not the physical width: for practical reasons the 
numerical computation of the Fourier transform in time in eq. (13) cannot integrate up to very 
large times 9 . Thus the width we see in the resulting plots is to a large extent contaminated by 
the fact that the time interval is finite in the numerical calculation (for the physical width to be 
visible unambiguously in these plots, the length of the time interval would have to be much larger 
than the lifetime of the quasi-particles). 

4 Occupation number 

4.1 Expression in terms of G + and G_+ 

Now that we know that at times x° > 100, the spectral content of the system reduces to simple 
quasi-particlcs, it makes sense to compute their occupation number. Recall that the creation and 
annihilation operators a k , ak are related to the field operator <j> via 



Bfe 

4 



i / d 3 x e ik ' x d x o 4>{x) 



= -i / d A x 



d x o 4>(x) . 



(34) 



9 One can check in the free case that this is a very singular Fourier transform. At k = 0, it is of the form 

J °° dt t exp(iwt) ~ <5'([^)- 
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From this, we get the following two reduction formulas 



a\a k ) 



ak.a k 



d 3 xd 3 y e *(»-i0 d x od y o G + „(x,y)\ x „ =y0 
d 3 xd 3 y e ik < x -ri d x od v o G-+{x,y)\ x0=y0 



(35) 



with the understanding that the times a; and y° are set equal only after the derivatives have been 
evaluated. 

It turns out to be more straightforward to calculate the sum of these two expectation values, 



i{a k + a k a k ^= d 3 xd 3 ye ik{x - y) l) x o*dyO G s {x,y)\ x0 __ 



(36) 



where G s = G^ + G |_, because in our framework the symmetric propagator G s is easier to 

compute than the separate G±^. The occupation number f k is related to the left hand side of 
eq. (36) by 

2w fc F(l + 2/ fe ) = (a{a k + a h a{) , (37) 

where V is the volume of the system and oj k the dispersion relation of the quasi-particles. 

4.2 Calculation of G s 

Let us now see how to compute the symmetric propagator G s (x,y) at LO, NLO and in the re- 
summation scheme we have developed to cure the pathologies related to secular divergences. At 
leading order, it is simply given by the product of two classical fields at the points x and y, 

G™(x,y) = 2ip(x)<p(y). (38) 

At next to leading order, G s is made of two pieces: 

i. a 1-loop correction j3 to one of the factors ip of the LO result, 

ii. a connected (tree-level) contribution Q s = Q^ + Q |_ that links the points x and y, 

G- LO (x, y) = 2 [ftxMv) + y{x)P{y)] + G s (x, y) . (39) 

The second term is given by 

d 3 k 



Qs{x,y) 



(27r) 3 2fc 



a +k (x)a- k {y) + a_ k (x)a +k (y) 



The a± k s can be formally related to the classical field tp(x) by (see eq. (5)) 

a±k(x) = / d 3 u [a±k ■ T u ] <p(x) , 



(40) 



(41) 



while for the tadpole /3 we can use eq. (21). Then, it is straightforward to combine the two terms 
to obtain 



G^°(x,y) 



d 3 u (3 ■ T u + i / d 3 ud 3 v f Jy 2k [a+k ■ TJ [a_ fc • T„] 



13 



From here, it is clear that one can perform the same resummation, where one exponentiates the 
quadratic part of the operator in the square brackets. This amounts to an average over Gaussian 
fluctuations of the initial classical field at a; = 0, 



G 



resummed/ 



(»,!/) = 2 / [DaDa]F[a,a] y{x)v{y) 



tpo+a 



(43) 



where the subscript ipo + a indicates the initial condition used at a; = to start the evolution of 
the classical field tp. F[a, a] is the Gaussian distribution of fluctuations defined in eqs. (7) and (8). 

4.3 Time evolution of /*. 

By combining the previous results, the occupation number obtained in this resummation scheme 
can be written as 



fk 



1 



1 



2uj k V 



[DaDa] F[a,a] 



d 3 x e lkx (ip(x°,x) + iuj k (p(x°, x)) 



(44) 



Vo+a 



In the evaluation of this formula, we use for the energy uj k = v k + m 2 with the mass fitted in 
the previous section (thus, we use a different mass at each time x°). The result of this calculation 
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Figure 5: Occupation number /& at various times in the evolution of the system. The grey band 
represents a fit by a Bosc-Einstcin distribution. The dashed red band is a fit by a pure power law 
k~ 5 ' 3 . The thin black line is a fit by a distribution of the form given in eq. (46). 

is displayed in the figure 5, where we show the occupation number at various stages of the time 
evolution, as well as three fits that we shall discuss shortly. 

Let us first briefly describe the main stages of the time evolution. At the initial time t = 0, 
only the zero mode is occupied and the higher modes have a negligible occupation number. This 
is a direct consequence of our setup, where the classical field is initially driven by a spatially 
homogeneous source. Then, shortly afterwards (this is already visible in the spectrum at t = 20) 
one sees an increase of the occupation in the non-zero modes, concomitant with a decrease of the 
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occupation in the zero mode (barely visible in the figure, due to the logarithmic vertical scale). 
The increase of the non-zero modes is most pronounced in a narrow band of k, where it peaks 
more than an order of magnitude above the rest of the curve. One can check 10 that this band of 
k coincides with the band of parametric resonance that we have discussed in detail in [48]. Thus, 
it appears that the dominant physics at early times is that of resonance, which leads to a quick 
increase of the occupation number in a narrow region of k. After t = 1000, the resonance peak has 
disappeared and the evolution becomes fairly slow. 

Let us now discuss fits of the occupation number, that are represented in the figure 5. The first 
two are a fit by a Bose-Einstein distribution, 



LW 



1 



op{^h—y) — l 



and a fit by a classical distribution of the form 



/class (,<V — 



T 



LU k - fl 



(45) 



(46) 



Interestingly, the best fit we could achieve with a Bose-Einstein distribution required a non-zero 
chemical potential. Although the particle number has no reason to be conserved in this theory 
(there is no symmetry protecting it), this suggests that changes of the particle number are slow 
compared to the evolution of the distribution in momentum space: a chemical potential at the latest 
times we have considered indicates that the particle number has not yet reached its equilibrium 
value (and its positive sign means that we have a particle excess). At weak coupling, this is rather 




Figure 6: Time evolution of the quasi-particle density in the system. Gray band: fit of the tail 
with a power law t^ 1 ' 4 . 

natural: inelastic processes have a much smaller rate than the elastic ones 11 , and therefore at 



10 See the appendix B of [48]. 

11 For the ifi 4 scalar theory that we consider here, a c i ~ c/ 4 , while o"i nc i ~ g 8 . The hierarchy between the elastic and 
inelastic time-scales is certainly less pronounced in QCD (sec [59]) and there it is unclear whether there is enough 
time for the formation of a transient state that has a non-zero chemical potential. 
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intermediate time scales the number of particles is an approximately conserved quantity. In order 
to check this hypothesis, we can evaluate the number density by summing the occupation number 
over all the modes k. This has been done in the figure 6. One sees indeed that, after a period of 
somewhat erratic evolution (that roughly corresponds to the time necessary to have well defined 
quasi-particles in the system), the number density decreases very slowly at late times, as a small 
negative power of time. 

It is also obvious from the figure 5 that a Bosc-Einstcin distribution does not fit well the 
occupation number in the tail at large k. In fact, the contrary would have been surprising, since our 
computation is essentially semi-classical. Naively, one may expect to obtain a classical distribution 
of the form T/(ujk — fx) (again, a non-zero /i is allowed if number changing processes are very slow), 
but one can check that such a distribution does not produce a better fit of the tail. At first sight, 
one could be tempted to blame this drop in the tail on the rarefaction of the lattice modes at large 
k (see for instance the figure 15 in [48]). However, this assumption does not hold if one does the 
same simulation with lower physical scales, as is done in the figure 7. On this figure, one sees the 
same drop, now occurring at a smaller value of k. If the drop was caused by lattice artifacts, one 
would expect it to occur at a fixed value of k (in lattice units) , no matter what the physical scales 
are. 
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Figure 7: Occupation number /& at various times in the evolution of the system, for a system 
initialized with a smaller energy density. 

It turns out that this drop has a rather trivial explanation. Firstly, note that a very good 
fit is obtained with the distribution given in eq. (46) (the thin black line in the figures 5 and 7), 
that differs from the naive classical distribution by an extra —1/2 term. This extra term, that 
makes the fit considerably better, has a simple origin: it comes from the —1/2 in the equation 
(44), which in the derivation of the formula for /& can be traced back to the non-zero commutator 
between creation and annihilation operators. Keeping this —1/2 correction in the definition of the 
occupation number for a semi-classical calculation is to a large extent an arbitrary choice. Indeed, 
such an approximation is expected to reproduce correctly the underlying quantum theory only 
in the region where the occupation number is sufficiently large. When this is the case, one has 
fk + 1/2 ~ fk, and therefore this 1/2 is not very significant. This also means that the drop of 
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the occupation number at k > 2 in the figure 5, although perfectly understandable in our semi- 
classical approximation, is obviously not a physical feature of the underlying quantum theory 12 . 
In our setup, there is one advantage in keeping the — 1/2 in eq. (44) though: if one does the same 
computation with a vanishing source J = 0, one gets identically /& = 0, which is of course the 
exact answer. Without this —1/2, one would have obtained /& = 1/2. 

4.4 Kolmogorov turbulence 

At t w 200, the modes in the resonance band reach their maximal occupancy, and start to subside 
afterwards, while the other non-zero modes continue to increase. While the resonance peak pro- 
gressively disappears, one sees in the figure 5 that the occupation curves tend to accumulate in the 
intermediate k range on a fixed line that is well fitted by a power law k~ 5 ' 3 . In this regime, the 
zero mode continues to decrease, while the occupation curve extends slowly into the hard region. 
Such a scaling with an exponent —5/3 in the power law is well know in the physics of turbulence 
(see the first part of [62] for instance). Typically, in Kolmogorov's turbulence, the energy cascades 
to the hard modes from a source localized in the soft sector, with an intermediate stationary 
distribution in between, that follows a power law k~ 5 ' 3 . In our case, the zero mode plays the role 
of this source, since it was initially the only occupied mode. In contrast to the usual setup in the 
study of Kolmogorov's turbulence [62], our system is closed and eventually the zero mode will run 
out of energy and will not be able to feed the cascade anymore. However, in our simulation, we 
have not reached the time at which this starts to happen. 

4.5 Bose-Einstein condensation 

In the figure 5, we saw that the occupation number at late times is best fitted by a distribution of 
the form of eq. (46). This fit however calls for two comments: 

i. the occupation number of the zero mode lies above the curve provided by this fit, 

ii. the best value of the chemical potential (p = 0.54) is close to the mass of the quasi-particles 
at this time, m = 0.58. 

These two seemingly unrelated facts have a common interpretation. As we have seen before, a 
chemical potential arises because the number of quasi-particles evolves very slowly in this system, 
and a positive /i is the reaction of the system to accommodate an excess of particles. However, it is 
clear from eqs. (45) and (46) that [i cannot be larger than the mass m - otherwise, the occupation 
number would become negative near k = 0. But having an upper bound on the chemical potential 
implies an upper bound on the particle density that these distributions can describe. What if 
the particle excess in the system is so large that the density is larger than this upper bound? 
When this happens, the excess of particles condenses on the zero mode, a phenomenon known as 
Bose-Einstein condensation. Dynamically, the system evolves towards a distribution made of two 



12 Interestingly however, the —1/2 term in eq. (46) is nothing but the second term in the expansion of the Bose- 
Einstein distribution in powers of (ui k — fJ.)/T. Therefore, at a formal level, keeping this 1/2 correction in the 
present semi-classical computation gives a better approximation of the full quantum theory. This point was already 
discussed extensively in [60, 61] in the context of the Boltzmann equation. 
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fk 



vf3(uJk- 



-1 



foS(k), 



(47) 



i.e. the chemical potential settles to the maximal value \i = m, and the extra particles go into the 
zero mode 14 k = 0. As one can see from the fit of the figure 5, the occupation number at low k 
appears to be precisely of the form of eq. (47). 

From the knowledge of the occupation number, it is easy to compute what fraction of the 
number of particles and what fraction of the total energy are contained in the zero mode at various 
times. This information is provided in the two plots of the figure 8. At early times (up to t ~ 100), 





Figure 8: Left: fraction of particles contained in the modes / < 
evolution. Right: fraction of energy contained in the modes \l\ < 



| fe | , at various stages of the time 
Ifcl. 



all the particles and all the energy are stored in the zero mode, as a consequence of our initial 
condition. At intermediate times (e.g. at t — 200), a large fraction of the energy is still in the 
zero mode, and the remainder is almost entirely in the resonance band. At the latest time we have 
considered (x° — 10 4 lattice units), the zero mode still contains about 35% of the particles and 
15% of the energy. 

On could argue that our computation does not demonstrate Bose-Einstein condensation, be- 
cause we started from an initial condition in which all the energy is already stored in the zero mode. 
What if we had started from a situation where the zero mode is empty? We have done that in the 
figure 9, in which the energy of the system is initially contained in the modes 15 (k x , k yi k z ) = (1, 1, 0) 
and (—1, —1,0) (the total energy being exactly the same as in the figure 5). Thus, at the initial 
time, the occupation number has a delta peak at a single, non-zero, energy. But then, one sees 
that a non-zero occupation number develops in the zero mode (and in other modes as well), to 
reach very large values in a rather short time. After this rapid transient regime, all trace of the 
original peak has been washed out. At late times, the distribution has become identical to the 
one encountered in the figure 5: all the modes except k = are described by a function of the 
form of eq. (46), and there is a particle excess in the zero mode. This study strongly suggests that 



13 Here we have written the quantum version of the distribution, but it has an analogue in our semi-classical 
approximation, where the first term is replaced by /dassC 8 )- This is discussed in more detail in the appendix A. 

14 Using the Boltzmann equation, it is easy to see that k = is the only mode where the extra particles can go. 
If the particles in excess occupy non-zero modes, then one does not have a fixed point of the Boltzmann equation. 
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We initialize the system in two opposite modes so that there is no net momentum in the system. 
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Figure 9: Occupation number /& at 
(k x ,k y ,k z ) = (1,1,0) and (-1,-1,0). 

times. 



various times, for a system initialized in the modes 
In the top right inset, we show the behavior at short 



Bose-Einstein condensation indeed occurs in this system 16 , when it is initially over-occupied. 

5 Further tests of the quasi-particle description 

5.1 Quasi-particle mass 

From the occupation number, we can further test the quasi-particle description of the system. 

A simple check is to compute the medium-induced mass of the quasi-particlcs, assuming that 

perturbation theory applies. The Hard Thermal Loop contribution to this mass has already been 

given in eq. (33), and we have represented this quantity in the figure 4 (blue curve), along with the 

mass obtained by fitting the location of the peak in the spectral function. In the region where quasi- 

particles are well defined (y° > 100), we see that the HTL value of the mass is systematically lower 

than the observed one. We can improve this result by including also the 1-loop vacuum contribution 

to the mass. At 1-loop in a <fi 4 theory, this is given by a tadpole graph whose expression can be 

written as 

d 3 k 



n-ioop 



{2ir) 3 2k 



(\ + h 



(48) 



Including the vacuum contribution to the mass improves significantly (see the red curve in the 
figure 4) the agreement between the theoretical prediction and the fit, indicating that higher-order 
corrections are presumably rather small for this value of the coupling (g = 1). Thus, it appears 
that the quasi-particle description is quite consistent: indeed, the occupation number computed 
from the fields themselves, when inserted into the 1-loop formula for the effective mass, reproduces 
very well the mass measured by fitting the peak in the spectral function. 



16 Of course, since the study is performed on a lattice, that has by definition a discrete spectrum, it is impossible 
to tell whether the distribution has a true S(k) term at the origin or whether it is a strongly peaked but otherwise 
regular function. 
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5.2 Residual interaction energy 

A quasi-particle description of a system is useful only if the residual interactions between the quasi- 
particles are weak - in other words, if the main effect of the interactions is simply to alter the 
properties of the particles (e.g. by generating an effective mass). This can be tested by computing 




Figure 10: Comparison between the actual energy density of the system (e) and the energy carried 
by the quasi-particles (e qp ) if one neglects their interactions. 

the energy density of the system by summing the energies of its quasi-particles, i.e. by assuming 
that they have no residual interactions , 



cqp 



d 3 k 

(2tt) 3 



fkVk 2 



(49) 



By comparing this quasi-particle energy to the actual energy density, e = (T 00 ), we can estimate 
the interaction energy of the quasi-particles and therefore the strength of their residual interactions. 
The result of this comparison is shown in the figure 10. We see that the true energy of the system is 
always below the energy of its quasi-particles, indicating that the residual interactions are attractive 
- which is indeed a standard result of a <fi 4 field theory. Moreover, as the time increases, the energy 
of the quasi-particles gets closer to the true energy, meaning that the quasi-particle description is 
better at late times. 



5.3 Entropy production 

From the occupation number, it is also possible to compute the entropy density, 

cPk 

(l + /fc)ln(l + /fc)-/ fc ln/ fc 



(2tt) 3 



(50) 



In our framework, replacing the true energy density by e qp is equivalent to substituting the expectation value 



of the interaction energy (V(ifi)) by \m 2 (<p 2 ), i.e. to a mean-field approximation 
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The time evolution of this quantity is shown in the figure 11 (green curve). One sees that the 
entropy density is multiplied roughly by a factor 20 during the evolution of the system (the initial 
value is low in our setup because the occupancy is entirely localized in the zero mode at t = 0). In 
the red curve, we have displayed the entropy that would have a gas of free bosons of equal energy 
density at thermal equilibrium 18 . The true entropy of the system gets close to the equilibrium 
entropy, but not exactly equal even at the largest times we have considered (the discrepancy 
remains of the order of 10 — 20%). This difference is presumably a combination of two factors: 
(i) the fact that the system is not yet fully equilibrated, and (ii) the residual interactions of its 
quasi-particles. 




Figure 11: Green curve (S): time evolution of the entropy density s defined in eq. (50). Red curve 
(S BE ): entropy of a non- interacting gas of bosons of equal energy density, in thermal equilibrium. 

The increase of the entropy defined by eq. (50) may seem paradoxical at first sight, given the way 
it has been obtained in our framework. Indeed, at the microscopic level, our system is described as 
an ensemble of classical field configurations that evolve according to the Euler-Lagrange equation 
of motion. This equation of motion is invariant under time reversal t — > —t, yet the entropy s 
computed in this system is clearly not invariant under this transformation - despite the fact that 
it is a functional of classical fields that have a time-reversible evolution. Let us recall here that 
the ensemble of classical field configurations that describe the system in our framework evolves 
according to the Liouville equation, 



d t F t + {Ft,H} = 0. 



(51) 



Instead of eq. (50), one could have defined an entropy based on the probability distribution Tt\ip, ip] 
of the field configurations in phase-space, 



S = — I [Dip Dip] Tt[p-,ft\ hi-Ft^, ip] 



(52) 



and it is easy to see that it is constant in time thanks to Liouville's theorem (see the appendix B). 
This means that, if one were able to determine the field configuration of the system at a given time, 



'Its small variations with time are due to the fact that the mass of the quasi-particles is not constant. 
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there would be no entropy increase because everything would be known about the microscopic state 
of the system. In contrast, the definition (50) is the appropriate definition when one knows only 
the single particle distribution in the system. Compared to cq. (52), a lot of information about 
the microscopic state of the system has been discarded by doing this. This coarse graining is the 
reason why the entropy given by (50) increases with time, while at the microscopic level the field 
configurations evolve via time reversible equations. 

6 Summary and outlook 

In this paper, we have pursued the study started in [48] of the time evolution of a system of scalar 
fields in a fixed volume box, in a rcsummation scheme that has been devised to cure the problem 
of secular divergences caused by parametric resonance. Although it is a rcsummation of quantum 
corrections that arise via loops, we have shown that this resummation can also be formulated as an 
average over classical field trajectories, with a Gaussian ensemble of initial conditions. This latter 
formulation is the one we adopt for a practical numerical implementation. In [48] , we have focused 
mostly on the evolution of the energy-momentum tensor of the system, and we have shown that 
the pressure relaxes towards its equilibrium value, thanks to a decoherence mechanism due to the 
fluctuating initial conditions. 

In the present paper, we have studied more microscopic properties of the system, in order to 
understand in more detail its time evolution. In particular, one question that was left unanswered 
in [48] is whether the system is in local thermal equilibrium at the time its pressure tensor reaches 
the equilibrium form. First of all, we have studied the existence of quasi-particles in the spectrum 
of the theory, by computing the spectral function of the system. Next, after having identified 
quasi-particle modes, we have computed their occupation number and the associated entropy. 

Let us first summarize the main results of this paper in a synthetic way, by displaying in parallel 
the time evolutions of the various quantities that we have considered separately so far. This is 
shown in the figure 12. From these plots, it appears that one can divide the time evolution in three 
stages that are qualitatively distinct 19 : 

i. < t < 100 : at these early times, the pressure of the system has not yet started to relax 
towards its equilibrium value p = e/3. Moreover, quasi-particles are not a good description 
of the system (their mass is not well defined, and there are extra spurious branches in the 
spectral function). The entropy displays only a very moderate growth during this era, the 
occupation number starts rising in the resonance band almost immediately, but the energy 
is still almost entirely contained in the zero mode, 

ii. 100 < t < 600 : this intermediate period starts when the occupancy in the resonance band 
reaches its maximum and starts to subside, while the occupation number rises in the other 
momentum modes. During this stage, the zero mode still contains the largest share of the 
total energy, and the remainder is contained predominantly in the resonant modes. The 
mains features of this era are the relaxation of the pressure towards its equilibrium value, 
and an important growth of the entropy. In this era, there are well defined quasi-particles, 
with a mass that is almost constant in time, 



19 The numerical values of the times quoted here are not absolute, but depend on the energy density in the system. 
Indeed, in a scale invariant theory, all the time-scales vary like e -1 ' 4 . Moreover, these time-scales depend on the 
coupling constant g 2 , and decrease as the coupling increases. 
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Figure 12: (Panels numbered 1 to 5 from the bottom to the top.) Panel 1: time evolution of 
the pressure. Panel 2: time evolution of the quasi-particle mass. Panel 3: time evolution of the 
entropy, compared to the entropy of a gas of same energy density at thermal equilibrium. Panel 
4: time evolution of the quasi-particle density. Panel 5: occupation number at various stages of 
the time evolution (the gray band is a fit at the latest time by a Bose-Einstein distribution with a 
chemical potential). 
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iii. 600 < t : in the late stages of the time evolution, the pressure is the equilibrium one, and the 
entropy displays only a marginal growth. However, the system is not yet fully equilibrated: 
the mass of the quasi-particles shows a clear decrease, while the occupation number continues 
to slowly expand towards higher momenta. During this stage, the occupation number shows 
some signs of Kolmogorov scaling, but perhaps even more compelling is the fact that it 
seems to evolve as would a system dominated by elastic collisions and which is ovcrpopulatcd 
compared to equilibrium - i.e. by developing a chemical potential equal to the mass and a 
Bose condensate at zero momentum. At the same time, the energy initially contained in the 
zero mode is progressively distributed among the higher modes. 

A very interesting observation is the appearance of a chemical potential, that we interpreted 
as resulting from a particle excess (compared to the value the particle density should have in 
equilibrium) combined to a fairly slow rate of inelastic processes. Eventually, the inelastic processes 
will wipe out the initial particle excess. However, because they are slow, there is an extended regime 
where the occupation number settles on a form that has a chemical potential. Moreover, for the 
bosonic statistics, the chemical potential cannot be larger than the mass of the quasi-particles. This 
implies that if the particle excess is too large, it cannot be accommodated solely by a chemical 
potential and the system stabilizes itself by the formation of a condensate at zero momentum. Our 
numerical results confirm this, and moreover indicate that this condensation occurs very early in 
the time evolution. Interestingly, the condition of overpopulation is also realized for the gluons 
produced initially in heavy ion collisions. Indeed, at a time Qj 1 (where Q s is the saturation 
momentum), the energy density is e ~ Qi/g 2 and the gluon number density is n ~ Ql/g 2 - Thus 
one has the dimensionless ratio ne -3 / 4 ~ g^ 1 ^ 2 ^S> 1, that should be of order 1 in chemical 
equilibrium. Depending on the strength of the number-changing processes, one may also expect 
the (transient) appearance of a gluonic chemical potential and the formation of a gluon condensate 
at zero momentum. 

One of the virtues of the formalism we have developed to address this problem is its seamless 
integration into the color glass condensate framework, since it is formulated in terms of classical 
field configurations (in fact, the resummation scheme we are using here was obtained as a by- 
product of the formalism developed for proving the initial state factorization of the large logs of 
l/.Ti.2). As a consequence, it can be adapted very easily to the Yang-Mills case, and to the study 
of thcrmalization in high energy heavy ion collisions. A crucial step in view of this application was 
performed in [47], where has been derived the spectrum of fluctuations at early times (i.e. before 
the Glasma instabilities start to develop) for gauge field fluctuations - in a system of coordinates 
and a gauge that are appropriate for numerical studies of the evolution of the glasma fields. 

Let us finally mention some possible connections with the idea of "eigenstate thermalization" 
proposed by Srcdnicki in [63] (see also [64, 65, 66]). This idea follows Berry's conjecture [67], 
which states that the high energy eigenstates of quantum systems whose classical counterpart is 
chaotic have extremely complicated wave-functions that for all practical purposes can be replaced 
by a random sum of plane waves (with Gaussian distributed coefficients). Assuming that this 
conjecture is satisfied, Srednicki proved that sufficiently inclusive measurements on such a state 
would lead to results that agree with thermal equilibrium. The main impediment to thermalization 
in a quantum system would thus be the fact that the system is usually not prepared in an energy 
eigenstate. For instance, for a system that is initially in a coherent state, thermal features would 
only become manifest once the various energy eigenstates superimposed in the coherent state have 
become incoherent. In the problem we have considered in the present paper and in [48], the initial 
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condition at t = appears to be closely related to the Wigner distribution of a coherent state and 
the associated classical dynamics is chaotic. Moreover, we have seen that the relaxation of the 
pressure is due to the loss of coherence between the various initial conditions. It would thus be 
very interesting to study whether there is a deeper connection between thcrmalization in quantum 
field theory and the ideas advocated by Srcdnicki. 
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A Evolution of the classical phase-space density 

In the figure 5, we have obtained a very good fit of the occupation number at late times by a 

function of the form 

T 1 

h = x • (53) 

Wfe - \x 2 

This fit works except for the zero mode, that is over occupied with respect to this distribution. 
We interpreted this distribution as the classical approximation of a Bose-Einstein distribution, and 
the —1/2 term was simply due to our definition of the occupation number. 

However, it is also interesting to forget the underlying quantum field theory we started from, and 
to consider in its own right the classical problem by which it is approximated. This reformulation 
is equivalent to solving the Liouvillc equation, 

d t T t + {T t ,H}=0, (54) 

given some Gaussian initial distribution. Therefore, if we adopt this point of view, we just have 
a (large) collection of coupled classical oscillators, and we follow their Hamiltonian flow in phase- 
space. In this appendix, we discuss some aspects of this classical dynamical system, that are 
relevant to the topics discussed in the rest of the paper. 

A.l Effective Hamiltonian 

Motivated by the observation of quasi-particles in the system, we may assume that there is a 
transformation of the fields and their conjugate momenta such that the Hamiltonian becomes a 
sum of quasi-free harmonic oscillators coupled only by weak residual interactions. In practice, this 
amounts to writing 



2V v r ' ) 4! 
d 3 x 1 L* + (V^) 2 + m V) + ^ - JmV ■ (55) 



2 V ' y r ' • r J • AV 2 
Wo U[ B 
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So far, we have just added and subtracted a mass term by hand, and the parameter m 2 is still 
arbitrary. In order to make the residual interactions small, one can choose the mean field value for 

2 

m , 

m 2 = ^-(<p 2 (x)) , ( 56 ) 

where the angle brackets denote an ensemble average 20 . The first part of this Hamiltonian can be 
rewritten as a sum of independent harmonic oscillators by going to Fourier space 21 , 

f d 3 k 1, ,2 1 2 | |2 

Ho= J (2^p # fc i ;gW ' (57) 

Ilk 

where cjfc = (k + m 2 ) 1 ' 2 and where tp k is the spatial Fourier transform of (p. 

This decomposition of the classical Hamiltonian into elementary harmonic oscillators plus resid- 
ual interactions is a good starting point to make connections with the study of the occupation 
number in the previous sections. Indeed, it is easy to check that eq. (44) is equivalent to 

- + / fe = ^. (58) 

In other words, (hk) is the occupancy of the mode k times u>k times the volume, plus a constant 
vacuum contribution Vu>k/2. This means that one should find a non-zero average value for (hk) 
even in the vacuum - actual particles in the mode k correspond to an excess of (hk) over VuJk/2. 

A. 2 Time evolution of the classical distribution 

In order to start the discussion, we first plot in the figure 13 the distribution J- t [tp, (p] of the classical 
field configurations at various stages of the evolution. Since we obviously cannot plot a functional, 
we have represented several Fourier slices: a slice being defined as the distribution of configurations 
in the plane 22 {\fwkjV\tpk\-, \v>k\/\ftJJkV)- In the figure 13, we have represented this distribution 
for three values of k: the zero mode, a mode in the resonance band, and a mode at some higher 
momentum. At the initial time, the zero mode of the fields is highly coherent, and its distribution 
is concentrated around a single point. In contrast, all the higher modes contain only fluctuations 
centered around the origin (0,0), with a width which is that of the vacuum fluctuations (i.e. the 
width that gives (hk) = Wfc/2). At the next time, t = 200, the zero mode has decohered and 
now fills almost uniformly a curve 23 of constant energy. We also see that the distribution of the 
resonant modes has expanded due to parametric resonance, while the harder modes still have the 
same distribution as the t = one. At later times, the expansion of the resonant modes reaches a 



20 One can check numerically that this mean field expression of the mass is in very good agreement with the 
measured mass of the quasi-particles. 

21 Note that this and subsequent integrals over d 3 k are strongly ultraviolet divergent. In all this appendix, one 
should have a lattice rcgularization in mind, so that the number of degrees of freedom in the system is finite. 

22 We rescale the Fourier components in this way so that the vacuum fluctuations look the same for all k. A value 
of (hh) proportional to lu^ corresponds to a circle of fixed radius (of order unity) in this plane. 

23 This is not exactly a curve of constant ho, as one can see from the fact that it is not a circle. This deviation 
from a circle is due to the fact that there is a large correction <ji 2 </Jq/4! coming from the interaction energy. Indeed, 
the zero mode, because of its large amplitude, is strongly interacting with itself. 
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Figure 13: Phase-space density in the (y/oJk/V\ipk\, \<fk\/V^kV) plane. From left to right: k = 0, 
1.07 (resonant mode), and 2. From top to bottom: t = 0, 200, 10 3 and 10 4 . Note the vastly 
different scale used for the zero mode (left column). The black circle represents the radius of the 
pure vacuum fluctuations. 
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maximum and then subsides, while the distribution of the hard modes begins to expand as well. 
Simultaneously, the zero mode distribution shrinks slowly, while remaining on a curve of constant 
energy (i.e. this seemingly constant energy is in fact slowly decreasing, due to a transfer of energy 
from the soft to the hard modes). 

To summarize these observations, the Fourier modes that initially have a large amplitude quickly 
decohere. If there are many such modes (instead of only one as in our numerical example), what 
happens then depends on whether the classical dynamics is chaotic or on the contrary intcgrablc: in 
a chaotic system, all these modes mix and fill uniformly an energy shell, while in the integrable case 
these modes would evolve independently and cover an invariant torus of much smaller dimension 
(for N modes, an energy shell is a manifold of dimension 2N — 1, while an invariant torus has 
dimension N only). On larger time scales, the energy carried by these initially large modes decreases 
slowly, and is transferred to the higher modes whose distribution expands as a consequence of this 
transfer. 

A. 3 Asymptotic behavior 

What is then the asymptotic distribution ^^[ipjip] of these classical fields? Let us assume first 
that the only quantity that is invariant under the Hamiltonian flow is the energy itself 24 . From 
the Liouville equation, it is clear that the asymptotic distribution must have a vanishing Poisson 
bracket with the Hamiltonian. This property is satisfied by any distribution that depends on 
ip, if only through the Hamiltonian H[tp, ip]. Such an asymptotic form for the distribution J^^ has 
strong implications on the average value of {hk). Indeed, if we can neglect the interaction energy 
(e.g. assuming that most of the interaction energy can be absorbed into an effective mass), then a 
distribution that depends only on H implies the equipartition of the energy among the modes, 

(hk) = const (independent of As) . (59) 

Note that this conclusion holds no matter what is the precise form of the function of W that J 7 ^ 

is equal to. When equipartition is in the classical phase-space is reached, the occupation number 

becomes of the form: 

const 1 

h = — r • (60) 

ujk I 

A. 4 Asymptotic behavior with constraints 

When we try to fit the late time curves of the figure 5 by an expression of the form (60), that 
lacks the chemical potential /n, the results are not good. This is best seen by plotting (hk) (see the 
figure 14), since eq. (60) would correspond to an horizontal line in this plot. In contrast, the much 
better fit of the figure 5 could be explained if the average value of hk at late times was of the form 

(hk) = const . (61) 

Wfc - \x 

This is obvious in the figure 14, where we obtain a good fit to (hk) by an expression of this form. 



24 The total spatial momentum of the system is also conserved, but it does not play any role here because we 
analyze the problem in the rest frame of the system. 
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Figure 14: Average value (hk) of the energy in the mode k. The brown band is a fit of the form 
Tujk/{ujk — /i), and its width reflects the expected statistical error in our Monte-Carlo simulation. 
The horizontal line shows the expected result if equipartition is reached. The red line represents 
the zero point energy Wfe/2. 

It is possible to understand eq. (61) if one assumes that the Hamiltonian flow not only conserves 
energy (exactly) , but also (to a good degree of approximation) the following quantity 25 



Af[tp, ip] 



d 3 k hk 



(2tt) 3 u k 



(62) 



It is obvious that this quantity has a vanishing Poisson bracket with the free quasi-particle part 
of the Hamiltonian TIq, and that only the residual quasi-particle interactions H[ nt can possibly 
make this quantity change. Thus, if the residual interactions are weak, this quantity should indeed 
be approximately conserved. If both % and Af are invariants, then the most general asymptotic 
solution of the Liouville equation must depend on ip, ip only through some combination of the form 
H — fiftf. For any distribution of this form, equipartition is replaced by 



hk — |U — ) = const (independent of fe) 



(63) 



which leads to eq. (61). Therefore, in the classical dynamical system, the quasi conservation of Af 
is what explains that we observed an energy distribution that differs considerably from the naive 
equipartition. 

Note that when k is large, eq. (63) is equivalent to (hk) = const. This implies that the constant 
in the right hand side must be positive. Considering now k — 0, this imposes /j, < m. Inserting 
now eq. (63) in the definition of Af, 



(Af) = 



d 3 k const 



(27r) 3 Uk — V 



(64) 



25 The average (Af) in the classical dynamical system is the counterpart of the number of quasi-particles in the 
quantum theory. 
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we see that it increases with /i, and reaches a finite (because the singularity at k = is integrable) 
maximum when /x = to. However, it may happen that the initial value of (M) is larger than this 
maximum. In this situation, the (/ifc)'s must be altered in order to accommodate this excess, but 
in such a way that eq. (63) remains valid. Any modification of (hk) for k ^ will violate eq. (63), 
so the only possibility is to modify (ho). Then, we see that the only way to change (ho) without 
violating eq. (63) is to have ijl = to. Thus, when there is an excess of (AA, the equilibrium value 
of (hk) takes the form 



(h k ) = A6(k) + B ■ 



^h 



u k 



(65) 



This explains why we found a chemical potential whose value is very close to the mass of the 
quasi-particlcs (within statistical errors). This phenomenon can be seen as an analogue in classical 
Hamiltonian dynamics of Bose condensation in quantum mechanics. 

A. 5 Do vacuum fluctuations thermalize? 

The fact that the Liouville evolution makes the distribution !F t evolve towards a classical thermal 
equilibrium leads to an interesting question: does the same happen if in our original problem there 
is no source coupled to the quantum fields, i.e. for pure vacuum fluctuations? In this case, the 
initial distribution is 

\<i>h\ 



Fo[<P,<p] =exp 



d 2 k \£k\ 2 + k 2 \<p k \ 2 



(2tt) 3 



(66) 



corresponding to (hk) = k/2. The consistency of our approach requires that the vacuum does 
not change over time. However, since eq. (66) is not a function of T-L, our previous considerations 
suggest that equipartition may also occur if we start from the pure vacuum fluctuations given by 
eq. (66). We have computed (hk) numerically for the pure vacuum case, by starting from the 



1 (pure vacuum) 




Figure 15: Average value of the energy hk in the mode k for pure vacuum fluctuations. The center 
of the brown band is u> k /2, and its width reflects the expected statistical error in our Monte-Carlo 
simulation). The horizontal line corresponds to equipartition. 



30 



distribution of eq. (66) and by evolving the corresponding classical field configurations also to a 
large time t = 10 4 . As shown in the figure 15, this computation shows no sign of equipartition of 
the vacuum energy, indicating that the vacuum is stable in our framework. 

This is at first sight an intriguing result. Indeed, how does the purely classical Liouville equation 
know that the distribution of eq. (66) represents the ground state of the corresponding quantum 
theory? In a sense, the Liouville equation seems to know more than what its "classical" qualifier 
might suggest. This result can be qualitatively explained by recalling a formal connection that 
exists between quantum mechanics and the classical Liouville equation. In quantum mechanics, a 
system can be described by a density operator p, whose evolution is driven by the Von Neumann 
equation, 

id t p+[p,H}=0. (67) 

It was noted by Wigner and Moyal that this quantum mechanical equation can be formulated 
equivalently as an evolution in the classical phase-space (see [68] for a recent review). To do this, 
one should first introduce the Wigner distribution associated to the density operator p. For a 
single degree of freedom, this is defined as 

W(q,p) = I dx e lpx (q + x/2\p\q - x/2) . (68) 

The Von Neumann equation can then be shown to be equivalent to the Moyal equation for W, 

d t W + {{W,H}} = , (69) 

where {{•,•}} is the Moyal bracket, obtained as the Wigner transform of the commutator. At 
this stage, one has a formulation of quantum mechanics that involves only quantities defined on 
the classical phase-space (the quantum mechanical aspects arc hidden in the fact that the Moyal 
bracket depends on h). The connection to the classical Liouville equation arises via the following 
two properties, 

i. the Moyal bracket {{•,%}} is equal to the Poisson bracket {•,%} if the Hamiltonian % is 
quadratic in the coordinates and momenta, 

ii. for any Hamiltonian, one has 

{{A,B}} = {A,B} + 0(h 2 ). (70) 

If the quantum system is in its ground state 0), its density operator po = 0/(0 is invariant under 
the Von Neumann equation. Equivalently, the corresponding Wigner distribution Wq is invariant 
under the Moyal equation. Since J-"o is the Wigner distribution of the same vacuum state, it is 
thus natural that it is left invariant by the Liouville evolution, since it is the h — > limit of the 
invariant Moyal evolution. 

B Liouville equation 

In this appendix, wc recall the derivation of Liouville's equation and some of its well known 
properties. 
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B.l Hamilton's equations 

In this appendix, we denote generically the canonical coordinates by Q, the corresponding canonical 
momenta by P, and the Hamiltonian by %. Hamilton's equations read 

Q=™ P = -™ (71) 



Let us introduce a few useful notations 



x,(% 



*-($?). ".(«</<£,). ™ 



Hamilton's equations can now be rewritten in the following compact form 

X = V . (73) 

This notation is suggestive of the fact that V is the velocity field induced on phase-space by the 
Hamiltonian. A crucial property of Hamiltonian flows is that they are incompressible, 

V ■ V = . (74) 

This identity is the essence of Liouville's equation. 

B.2 Liouville's equation 

Consider an ensemble of such dynamical systems, all described by the same Hamiltonian %. At 
the time £, their distribution in phase-space is described by a density J- t [Q,P], and we wish to 
derive an equation that describes the time evolution of this distribution. Naturally, the number 
of systems in the ensemble is not changing since each system evolves independently of the others. 
Thus, Ft is the density for a locally conserved quantity. We have seen before that each point 
in phase-space moves with the velocity V. Therefore, we can write a continuity equation that 
expresses this conservation, 

d t T t + V-(T t V) = 0. (75) 

By using eq. (74), we can rewrite this equation as 26 

d t + V-V^ t = 0. (76) 

Note that the term V ■ VJ-j is nothing but the Poisson bracket {Ft,T-L}- Thus, we see that both 
eqs (75) and (76) are equivalent to the usual form of the Liouville equation, 

d t T t + {F t ,H} = 0. (77) 

These alternate forms of the Liouville equation are very useful, because they recall us its origin 
as the continuity equation for the density of systems in phase-space, and because they make some 
properties of Hamiltonian flows more transparent - in particular thanks to the Row derivative D t 
introduced in eq. (76). 



26 In the form DtTt = 0, the Liouville equation is equivalent to the Liouville theorem, that states that Tt is 
constant along the Hamiltonian flow lines. 
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B.3 Basic properties 

In this section, we derive some elementary properties of Liouville's equation. To state some of 
these properties, it will be useful to denote [dT] the measure on phase-space. 
The most elementary property, 

/ [dT] F t = const , (78) 

is in fact just the integral version of Liouville's equation itself. It is simply another way of stating 
that the number of systems in the ensemble does not change over time. Note that since Ft is 
constant along the flow lines, the same is true of any local function F(F t ). Thus, eq. (78) can be 
generalized by replacing under the integral Ft by any function F(Ft). A similar statement can be 
made about the energy if we note that 27 

D t H = . (79) 

(which means that the Hamiltonian does not vary if we follow the flow). Then, by multiplying 
eq. (76) by H, and by integrating over phase-space, we get 

[[dT] F t H = const . (80) 

This equation simply says that the total energy of our ensemble of systems is conserved. 
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